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ABSTRACT 

In this second paper in our series, we continue to test primordial scenarios of globular cluster formation 
which predict that globular clusters formed in the early universe in the potential of dark matter minihalos. In 
this paper we use high-resolution A^-body simulations to model tidal stripping experienced by primordial dark- 
matter dominated globular clusters in the static gravitational potential of the host dwarf galaxy. We test both 
cuspy Navarro-Frenk- White (NFW) and flat-core Burkert models of dark matter halos. Our primordial globular 
cluster with an NFW dark matter halo survives severe tidal stripping, and after 10 orbits is still dominated by 
dark matter in its outskirts. Our cluster with Burkert dark matter halo loses almost all its dark matter to tidal 
stripping, and starts losing stars at the end of our simulations. The results of this paper reinforce our conclusion 
in Paper I that current observations of globular clusters are consistent with the primordial picture of globular 
cluster formation. 

Subject headings: globular clusters: general — methods: N-body simulations — dark matter — early universe 



1. INTRODUCTION 

The primordial scenario of glob ular cluster (GC ) forma- 
tion was first proposed by Peebles & Dicke] Jl968h . Their 
paper considers GCs forming in the early universe out of 
primordial gas fluctuations on a Jeans mass scale. Realiza- 
tion that the dominant mass component of the univ erse is 
dark m atter (DM) led to the revision of this idea in iPeeblesI 
(QUI, where GCs are proposed to have formed in the po- 
tential wells of DM minihalos in the early universe. In 
this revised picture GCs can be considered as some of the 
earliest galaxies in the universe, with a baryons dominated 
core and extended DM halo. This picture of DM-dominated 
GCs forming in the earl y universe was later consi dered 
by ma ny authors, includin g Rosenblatt, Faber, & Blu menthall 
(119881). IFadoan Jimenez. & Jonesl d!997h ICenl (12001). 
Bro mm & Clarke! <2002l> . and lBeaslev et alJ ll2003l) . 

Primordial GCs could have formed somewhere between a 
redshift of ~ 30, when the first stars in the universe are be- 
lieved to have b een born inside m inihalos with total masses 
10 s - 10 6 M Q < Yoshi da et alJ2003l) . and redshift of ~ 6, when 
the re ionization of the universe was complete ([Becker et alJ 
120011) . It is often assumed that the very first stellar objects 
produced enough Lyman- Werner photons (capable of dissolv- 
ing H2 molecules, which are the main coolant in the pris- 
tine gas) to prohibit star formation in most of the low mass 
(< 10 8 Mq) gas-ric h halos. Cosmological simulati ons with 
radiative transfer of Ricotti. Gnedin. & Shul| ll2002ft suggest 
another possibility: in their model, numerous mini-galaxies 
with masses of < 10 8 Mq form inside relic cosmological H II 
regions. The observed positive feedback can be explained by 
high non-equilibrium fraction of free electrons in defunct H II 
regions, which leads to enhanced H2 production and hence 
star formation. This is an interesting mechanism for forming 
GCs, as it naturally explains the observed higher specific fre- 
quency of GCs in higher density regions of the universe (such 
as around cD galaxies at the center of clusters of galaxies), 
where large cosmological H II regions should have formed in 
the early universe. 

In a widely accepted hierarchical paradigm of structure for- 
mation in the universe, small bound objects form first, later 
merging to form objects of increasingly larger size. In this 



picture, primordial GCs would first experience a major merger 
with comparable mass halos, or be accreted by a larger halo 
of a dwarf galaxy size. In their simulations of the forma- 
tion of a dwarf g alaxy in the early universe (at z > 10), 
Bromm & Clarke (2002) observed a formation of GC-type 
baryonic objects inside DM minihalos, which later merged 
to form a dwarf galaxy. By the end of the simulations, bary- 
onic cores, believed to repres ent proto-GCs, had appa rently 
lost their individual DM halos. iBromm & Clarkd ( 120021) sug- 
gested that DM was lost because of the violent relaxation ac- 
companying the major merger. The lack of resolution did not 
allow the authors to reach more definitive conclusion as to the 
fate of DM in their proto-GCs. 

In the next step, dwarf galaxies containing primordial GCs 
would merge to produce objects of increasingly larger size. 
This process should continue at the present time. A well- 
known example is that of Sagittarius dwarf galaxy, which is 
presently in the process of being disrupted by the tidal forces 
of Milky Way and depositing i ts system of GCs into the halo 
of our Galaxy (iBellazzini. Ferraro. & Ibatal20 03). 

In the primordial picture of GC formation, GCs can 
be considered as small nucleated dwarfs accreted at 
some point by larger galaxies, which lose their extra- 
nuclear material (baryons and DM) due to stripping 
by the tidal forces of the accreting galaxy. Such 
mechanism was invoked to explai n the formation of 
a few of the larg est GCs: M5 1 (lLavden & Saraiedinil 
2000), a) Centaur jjjTsjtdriva^in escu. & K orchagin ]2003t 
Mizutani. C hiba. & Sakamoto 20031). and Gl (iBekki & Chibal 
2004), which is the satellite of the M31 galaxy. In the quoted 
papers, these GCs are considered to be exceptions. We sug- 
gest instead that this method of forming these GCs is a rule 
rather than exception, and is fully consistent with primordial 
GC formation scenarios. In this picture, the above largest 
GCs are examples of the most recent GC-producing accretion 
events, and represent a small fraction of the total number of 
proto-GCs accreted directly by a large galaxy (Milky Way or 
M31) in the modern universe (thus skipping the intermediate 
step of being accreted by a larger dwarf galaxy in the early 
u niverse) . 

iMoord ( Ti99f3) used numerical simulations of a DM- 
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dominated GC orbiting inside the static potential of Milky 
Way halo and ob servations of tidal tails of M2 of 
iGrillmair et alJ (I 19951) to claim that primordial scenarios of 
GC formation are inconsiste nt with observations of G alac- 
tic GCs. As discu ssed by iBromm & Clarke! J2002I) and 
Beaslev_£LaiJ{2003), the model of Moore was too simplistic 
to address this issue. Most importantly, the author considered 
all Galactic GCs being directly accreted by the Milky Way, 
which is definitely not the case in the hierarchical picture of 
structure format ion in the un iverse. Moreover, the DM sub- 
halo model of Moore ( 1996) is a lowered isothermal sphere 
with arbitrary (not cosmologically normalized) structural pa- 
rameters, which makes it very different from the DM halos 
expected to host a GC in the early universe. 

One also has to be careful with claims on the existence 
of "extratidal features" around Galactic satellites, which in- 
cludes both GCs and dwarf sp heroidal galaxies. In the case of 
u Centauri. lLaw et a l](|2003b showed that the appare nt "tidal 
tails" observed by Leon, Mevlan, & Combes (2000) around 
this GC are most probably caused by spatially variable fore- 
ground dust extinction in its outskirts. In the case of the Draco 
dwarf spheroidal galaxy, the detecti on of apparent "extrati- 
dal fe atures" around this galaxy by Irwin & Hatzidimitriou 
( 1995) was later refuted by observations of Odenkirchen et al. 
(2001) which had much higher sensitivity. The "extratidal" 
extension of the radial surface brightness profiles of Galactic 
GCs beyond a "tidal" King radius is almost always seen at 
the levels lower than the inferred background contamination 
level, and could be an artifact of the background correction 
procedure. One beautiful e xample of a GC w ith obvious tidal 
tails is that of Palomar 5 ( Odenkirchen et al. 2003). As we 
will show in this paper, a presence of tidal tails in some GCs 
is not inconsistent with these GCs having formed with DM 
halos in the early universe, and hence with primordial scenar- 
ios of GC formation. 

In this series of two papers, we are using high-resolution 
Af-body simulations of GCs with DM halos to test primor- 
dial scenarios of GC formation. In particular, we are trying 
to answer the following questions: 1) Are there obvious sig- 
natures of DM presence in GCs? 2) Will DM in a primor- 
dial GC survive hierarchical process of assembling galaxies 
in the e arly universe? In the first paper, Mashchenko & Sills 
(2004b, hereafter Paper I), we considered the initial relax- 
ation of a stellar cluster at the center of a DM minihalo in 
the early universe (around z = 7). The structural parameters 
of DM halos were fixed by results of cosmological ACDM 
simulations. The initia l non-equilibrium con figuration of stel- 
lar clusters was that of Mashchenk o & SiirsI fc004al hereafter 
MS04), where we showed that purely stellar (no DM) ho- 
mogeneous isothermal spheres with the universal values of 
the initial density p,- * = 14 Mq pc~ 3 and velocity dispersion 
er,.» = 1.91 km s -1 collapse to form GC-like clusters, with 
all bivariate correlations between structural and dynamic pa- 
rameters of Galactic GCs being accurately reproduced. We 
showed in Paper I that many observational features of Galactic 
GCs, used traditionally to argue that these systems are tidally 
limited purely stellar clusters without DM, can be produced 
by the presence of significant amounts of DM in their out- 
skirts. In particular, in warm collapse models (with the initial 
virial parameter for the stellar core of 0.27 < v < 1.7) we 
observed the formation of an apparent "tidal" cutoff in sur- 
face brightness profiles, with no unusual features in the outer 
parts of the velocity dispersion profile. In cold collapse mod- 



els (with v < 0.27), on the contrary, an apparent "break" in 
the outer parts of both surface brightness and velocity disper- 
sion profiles is seen, which can be mistakingly interpreted as 
presence of "extratidal" stars heated by the tidal field of the 
host galaxy. 

The results of Paper I are directly applicable to dynamically 
young intergalactic GCs and GCs which have not experienced 
significant tidal stripping in the potential of the host galaxies. 
In this second paper, we test the regime of severe tidal strip- 
ping of our primordial DM-dominated GCs in the potential of 
the host dwarf galaxy. As we discussed above, being accreted 
by a dwarf galaxy would be a typical fate of primordial GCs 
formed at high redshift. We use warm collapse models of pri- 
mordial GCs from Paper I to set up the initial conditions for 
the present paper. We c o nsider both simulations suggested 
Navarro. Frenk. & White ( 1997, hereafter NFW) and obser- 
vationallv motivated IB urkertl (119951) density profiles for our 
DM halos. 

This paper is organized as follows. In Section |2] we de- 
scribe our models. In Section[3]we discuss the results of the 
simulations. We conclude with Section|4] 

2. MODEL 
2. 1 . Physical Parameters of the Models 

In our model, we explore the impact of tidal fields on DM- 
dominated GCs orbiting inside a static potential of host dwarf 
galaxy with a virial mass of 10 y Mq, GCs are assumed to be 
accreted by the dwarf galaxy at a redshift of z = 3. Structural 
parameters (virial radius R^, scale radius R s , and concentra- 
tion C = R v i r /R s ) of the host galaxy are fixed by the results 
of cosmological ACDM simulations, and at z = 3 are equal to 
R vil = 8.16 kpc, R s = 1.21 kpc, and C = 6.75 (see eqs. [1-9] 
from Paper I). Throughout this paper we assume the follow- 
ing values of the cosmologic al parameters : fi m = 0.27 and 
H = 11 km s" 1 Mpc" 1 (Sper eel etall l2003). Galaxies with 
M v ; r = 10 9 Mq correspond to ~ la de nsity fluctuations col- 
lapsing at z = 3 (Jgarkana & Loebl20 01. their Fig. 6) and hence 
are very common at that redshift. 

The initial parameters of our GCs are identical to those of 
our warm collapse models from Paper I. We could not use the 
cold collapse model of Paper I because it would be an impos- 
sible task with the present day technology with the method we 
are using (A^-body simulations). Indeed, all the models pre- 
sented in the present paper for the warm collapse only case 
required 3.6 CPU-years to run. The cold collapse models 
from Paper I (models C„j,) took ~ 7 times longer to run than 
the warm collapse ones (models W Wj /,) for the same evolution 
time, which is mainly due to ~ 6 times shorter stellar crossing 
time in the cold models. In addition, to run the cold collapse 
models of Paper I for the much longer evolution time of the 
present paper (3 Gyr), we would have to increase significantly 
the accuracy of the integration, and increase the number of 
DM particles by a factor of 20 (up to 10 7 particles) to avoid 
artificial mass segregation effects which should be very sig- 
nificant for such long evolution time. As a result, the time 
required to simulate the cold collapse models for 3 Gyr would 
be two or more orders of magnitude longer than for the warm 
collapse case, or a few hundred CPU-years, which is obvi- 
ously not feasible. The reasons for choosing the warm versus 
hot collapse model are that the warm one corresponds to a 
more typical GC (with the mass of ~ 10 5 M Q for the warm 
model versus ~ 10 4 M Q for the hot one), and that the warm 
model of Paper I exhibits interesting tidal-like features in the 
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outskirts of the cluster (namely, a cutoff in the surface bright- 
ness profile) caused by the presence of DM halo, whereas the 
hot collapse model does not have any such features. 

An isothermal homogeneous stellar sphere with the univer- 
sal values of the initial stellar density and velocity dispersion, 
p,> = 14 M Q pc" 3 and cr (> = 1.91 km s" 1 (MS04), is set at 
the center of a DM halo with either NFW or Burkert density 
profiles. The virial mass of the DM halo is mdm = 10 7 M Q . 
The stellar mass is m* = 8.8 x 10 4 M Q , so the baryonic- 
to-DM mass ratio is x = 0.0088. Our fiducial value of x, 
on one hand, is larger than the fraction of baryons in GCs 
of ~ 0.0025 jMcLaug hlinl[T999i. and on the other hand, is 
smaller than the universal baryonic-to-DM density ratio of 
~ 0.20 ( Snergel et alJl2003l) . We assumed GCs have formed 
at z = 7, so the structural parameters of their DM halos are 
r vir = 885 pc, r s = 181 pc, and c = 4.88 (from eqs. [1-9] of 
Paper I). 

The initial virial ratio for our stellar clusters is z/» = 0.54. 
The initial stellar radius is 11.2 pc. As we showed in Pa- 
per I, in the absence of DM our stellar cluster will first 
collapse (with the smallest value of the half-mass radius of 
'"♦.min = 3.0 pc), and then bounce to form a relaxed cluster with 
a flat core, with the radius of ro = 3.0 pc, and a steep power- 
law outer density profile. The crossing time at the half-mass 
radius of the relaxed cluster is t* = 0.49 Myr. The crossing 
times at the half-mass radius for our DM halos are 52 and 
61 Myr for NFW and Burkert profiles, respectively. 

For consistency, GCs with NFW halos are assumed to orbit 
a host galaxy with NFW profile, and GCs with Burkert halos 
orbit a Burkert host galaxy. In all models, the pericentric dis- 
tance is R p = R s /2 = 0.60 kpc, and the apocentric distance is 
Ra = Rs x 2.5 = 3.02 kpc. We chose such a small value of R p 
to explore a regime of strong tidal stripping, potentially ap- 
proaching a violent relaxatio n mode of DM str ipping in the 
GC formation simulations of Bromm & Clarke! ( 120021) . Our 
apocentric-to-pericentric distance ratio R a /R p = 5 is similar 
to orbits of substructure in cosmological ACDM simulations. 
The theoretical radial orbital period is P = 320 Myr for the 
host galaxy with NFW profile, and P = 340 Myr for the Burk- 
ert case. In our simulations, subhalos experience dynamic 
friction caused by the motion of the remnant in the halo of 
tidally stripped DM particles. This effect is the most notice- 
able during first 1-2 orbits (for the following orbits the mass 
of the subhalo becomes too small to be affected by the dy- 
namic friction), and leads to slightly smaller values of R p , R a , 
and P. All our simulations are run for t2 = 3 Gyr. (This cor- 
responds to the range of redshifts of z = 3 ... 1.2.) The actual 
radial orbital periods in our models are P = 270 Myr for NFW 
potential and 290 Myr for Burkert potential, so the total num- 
ber of orbits is 10-11. 

2.2. Numerical Parameters of the Models 

In total, we run 9 different simulations (see Table[Q. Mod- 
els S (stars only), D„j, (DM only), and SD„ (DM + stars) are 
evolved in isolation (no external static gravitational field), and 
are used as reference cases for the models DO„.b (DM only) 
and SDO„.fc (DM + stars) where subhalos orbit the static po- 
tential of the host galaxy. (Here subscripts "n" and "b" denote 
NFW and Burkert DM profiles, respectively.) 

For models containing DM, we generate initial distribu- 
tion of DM particles using the rejection method described 
in Paper I. Halos are represented by 10 6 particles with in- 
dividu al masses of 10 and a softening leng th of €dm = 
1.5 pc. Kazantzidis, Mag orrian. & Moord ( 120 04) showed that 



TABLE 1 

Numerical parameters of the models 



Model 


N, 


6* 






t2 


Aw 


Note 






pc 




pc 


Gyr 


Gyr 




S 


I0 4 


0.30 






3 


2 X 1Q- 4 


1 


D„ 






10 6 


1.5 


3 


2 X 10~ 4 


2 


D/, 






LO 6 


1.5 


3 


2 X 10~ 4 


2 


SD„ 


10 4 


0.30 


10 6 


1.5 


3 


2 X 10^ 


3 


SD 6 


1() 4 


0.30 


10 6 


1.5 


3 


2 X 10~ 4 


3 


DO,, 






10 6 


1.5 


3 


2 X 10~ 4 


4 


D0 6 






1() 6 


1.5 


3 


2 x 10~ 4 


4 


SDO„ 


10 4 


0.30 


10 6 


1.5 


3 


2 x 10- 4 


5 


SDO fc 


10 4 


0.30 


10 6 


1.5 


3 


2 x 10- 4 


5 



NOTE. — 1) Stars only. 2) DM only. 3) Stars + DM halo. 4) DM only 
halos on orbit inside a static potential. 5) DM halos with stellar cores on orbit 
inside a static potential. Here /V* and Wdm are the number of stellar and DM 
particles; e* and edm are the softening lengths for stars and DM; t2 and A? max 
are the total evolution time and the maximum value for individual timesteps. 



for studies dealing with tidal disruption of substructure one 
should not use so-called "local Maxwellian approximation", 
where local velocity distribution is assumed to be multivariate 
Gaussian, to set up initial conditions, as such configuration is 
not in equilibrium for cuspy DM density profiles and can lead 
to artificially high disruption rate for subhalos. For our mod- 
els, we explicitly use phase-space distribution functions (DFs) 
to assign the components of velocity vectors to different par- 
ticles, which guarantees the central part of the halo to be in 
equilibrium initially. To calculate DFs, for NFW p rofile we 
use the analytical fitting formula of Widrow (2000), and for 
Burkert profile we use our own fitting formula (Paper I, Ap- 
pendix A). 

Our halos are truncated at a finite radius r y „, which results 
in the outer part s of the halos being not in equilibrium. As 
we will see in § 13.21 this effect has negligible impact on the 
results of our simulations. 

For models containing a stellar core, stars are repre- 
sented by 10 4 equal mass particles, with individual masses 
of 8.8 Mq. The softening length for stars is e = 0.30 pc. 
Stars are set up as a homogeneous isothermal sphere located 
at the center of the DM halo (for models containing DM). 
Stellar particles initially have a Maxwellian distribution of 
velocity vectors. As we showed in MS04, such initial non- 
equilibrium configuration of a stellar cluster leads to forma- 
tion of core-halo structure, with a radial density profile re- 
sembling that of GCs, after the initial relaxation phase. The 
model of MS04 also successfully reproduces all empiric bi- 
variate correlations between structural and dynamic parame- 
ters of Galactic GCs, given that all proto-GCs start with the 
same values of the stellar density p, > = 14 M Q pc" 3 and ve- 
locity dispersion cr,.* = 1.91 km s" 1 . We use these values of 
p^* and cr,> for setting up our models. 

Our models SDO„ and SDO;, are first evolved in isolation 
(with no static gravitational field) for 120 and 170 Myr, re- 
spectively. This allowed stars and DM at the center of the 
halo to reach a state of equilibrium. 

We use a parallel version of the mu ltistepping tree code 
GADGET (Springel . Yoshida & Whitell2001l) to run our sim- 
ulations. The values of the code parameters which control 
the accuracy of simulations are the same as for our warm 
collapse models W„ j from Paper I. In particular, we use a 
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very conservative (small) value of the parameter rj = 0.0025 
controlling the values of the individual timesteps, which are 
equal to (2rje/ 'a) 1 ? 2 , where a is the acceleration of a particle. 
Also, the individual timesteps are not allowed to be larger than 
Af max = 2 x 10~ 4 Gyr (see Table[Q. As a result, we achieved 
an acceptable level of accuracy in our simulations, with the 
total energy change AEtot being equal to 1.8% after 3 Gyr 
(or more than 6000 crossing times at the half-mass radius) for 
our purely stellar model S. For models with DM (D„^ and 
SD„ J,) this number is significantly smaller (AE tot < 0.07%), 
because their total energy budget is dominated by low den- 
sity parts of the DM halo. We cannot estimate AE tot for our 
static field models DO,,.;, and SDO,,.;,, but we expect the inte- 
gration accuracy for these models to be comparable with that 
of the corresponding models without the static field (D„ b and 
SD„, fc ). 

For all models, we output 100- 1000 snapshots for different 
moments of time. In every snapshot we identify a gravitation- 
ally bound structure (if present) using the same routine as in 
Paper I. This procedure consists of two main steps. 1) We 
use the program addqr avity, which is based on the algo- 
rithm o f lDehnenl 12000) and is part of the NEMO 1 software 
package, to assign local gravitational potential $ values to in- 
dividual particles (both DM and stars). We use the softening 
length values from Table ^for purely stellar and purely DM 
models, and an intermediate value of e = 1 pc for hybrid mod- 
els. Next, we find the 1000 (100 for model S) particles with 
lowest potential. For these particles we calculate weighted 
six phase-space components of the halo center, using a nor- 
malized potential ipt = (<J> max - $,)/($ ma x - $min) as a weight. 
Then we recenter the snapshot both spatially and in velocity 
to the halo center. 2) We remove all unbound particles (with 
velocity module v > [-2G&] '/ 2 ) in a single step and recom- 
pute the potential for the remnant using addgravity. We 
repeat this unbinding procedure in an iterative manner until 
there are no unbound particles. 

The above procedure is reasonably fast and accurate. It 
failed to find a bound subhalo only in the second half of DO„ 
run. For a few "failed" snapshots we had to use the program 
SKID 2 , which is significantly slower than our simple unbind- 
ing procedure. (SKID is often used to find gravitationally 
bound structure in the results of cosmological simulations.) 

3. RESULTS OF SIMULATIONS 

3.1. Long-Term Dynamic Evolution of the Stellar Cluster 

In the first two lines of Table[2]we show the model param- 
eters for an isolated stellar cluster without DM (model S) for 
two moments of time — soon after the initial relaxation of the 
cluster (t ~ 0.28 Gyr; first line) and at the end of the simula- 
tions (t ~ 2.84 Gyr; second line). As you can see, some stellar 
cluster parameters (such as half-mass radius r^* and central 
velocity dispersion a c ) stay virtually the same throughout the 
simulations, whereas others (central density p c , central sur- 
face brightness £y, King core radius ro, and core mass frac- 
tion /o) undergo significant changes. 

On a more detailed level, these changes can be followed in 
Figure ^ where we show averaged radial density profiles for 
model S for four different time intervals. As you can see in 
this figure, the flat core of the cluster becomes smaller and 
denser with time. Interestingly, the "dent" feature, seen out- 
side of the core in the radial density profile of a freshly relaxed 
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FIG. 1 . — Averaged radial stellar density profiles for model S for four time 
intervals: t = 0. 12 . . . 0.37, 0.37 ... 0.87, 0.87 ... 1.61, and 1 .61 ... 3 Gyr (from 
bottom to top in the left part of the figure, and from left to right in the right 
part). 



cluster, disappears in significantly evolved clusters. The outer 
part of the profile becomes increasingly more shallow with 
time, with the following values of the power-law exponent: 
7 = -7.4, -5.9, -5.1, and -4.6. 

As you can see in Table [2] by the end of the simulations, 
~ 50% of core stars in model S evaporate. These stars stay 
gravitationally bound to the cluster, populating its outer parts 
and making the outer density profile more shallow. 

The observed secular evolution of our model stellar cluster 
is very similar to gravothermal instability (or core collapse) 
known to occur in real GCs. The code we use to simulate 
the cluster is collisionless: it softens gravitational potential of 
particles separated by less than the softening length e* (which 
is comparable to the average distance between particles in the 
cluster). As a result, close encounters between particles are 
not treated correctly. Moreover, our stellar particles do not 
represent individual stars, as they are ~ 10 times more mas- 
sive than stars in GCs. Nevertheless, it appears that our colli- 
sionless simulations capture the essence of the long-term dy- 
namic evolution of GCs, probably because the main driving 
mechanism for such evolution is not infrequent very close en- 
counters between stars, but rather the cum ulative effect of nu- 
merous weak interactions (Spitzer 1987), which are resolved 
reasonably well by our code. 

The similarity between our stellar model's secular evolu- 
tion and gravothermal instability is also seen on a more quan- 
titative level. The idealized model of core collapse of a GC 
predicts a correlation betwee n the central d ensity and the ra- 
dius of the core: p c cx t"q 221 lSpitzerlll987l) . In our model S, 
we observe a very similar correlation: p c cx r^ 20 for t he time 
interva l t = 0.3 ... 3 Gyr. Both the idealized theory of Spitzer 
( 1987) and our model exhibit very mild evolution of central 
velocity dispersion: a c cx and a c cx Tq 01 , respectively. 

The fact that our stellar particles are ~ 10 times more mas- 
sive than stars in GCs results in the pace of secular evolution 
in our models being significantly faster than in real GCs. To 
estimate how dynamically old our models are at the end of 
simulations we use the analytical theory of core collapse of 
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TABLE 2 

Model parameters at the end of the simulations 



Model 


rh* 




Pc 


Rh,* 


Rhb 


CO 




ro fo 


T 


"'DM 


*Vi.DM 


r p 


r m 




Mq pc 


km s _1 


Mq pc 3 


pc 


pc 


km s - ' 


mag arcsec -2 


pc 


MqL- q 1 


Mq 


pc 


pc 


pc 


S a 


8.80 x 10 4 4.73 


4.02 


360 


3.65 


2.55 


3.82 


18.66 


2.75 0.211 


1.47 










s 


8.78 X 10 4 4.68 


4.05 


1300 


3.56 


1.33 


3.96 


17.91 


1.43 0.105 


1.53 










D„ 




















9.68 x 10 6 


434 






D/> 




















9.66 x 10 6 


497 






SD„ 


8.80 X 10 4 4.31 


4.46 


980 


3.26 


1.60 


4.41 


18.05 


1.83 0.153 


1.79 


9.70 x 10 6 


424 


7.7 


19.1 


SD b 


8.80 x 10 4 4.77 


4.05 


1400 


3.59 


1.31 


3.98 


17.91 


1.42 0.105 


1.56 


9.68 x 10 6 


484 


17.8 


52.9 


DO,, 




















3.83 x 10 4 


28.4 






SDO„ 


8.80 x 10 4 4.31 


4.41 


1200 


3.27 


1.44 


4.37 


17.97 


1.67 0.138 


1.81 


2.67 x 10 5 


38.0 


8.6 


24.5 


SDO,, 


8.77 X 10 4 4.72 


4.01 


1400 


3.57 


1.28 


3.98 


17.89 


1.40 0.106 


1.56 


4.70 x 10 4 


65.1 


25.7 





NOTE. — Here m„, r/, cr c , and p c are mass, half-mass radius, central velocity dispersion, and central density of bound stellar clusters; Ri, „, <jq, and Ey 
are projected half-mass radius, half-brightness radius (where surface brightness is equal to one half of the central surface brightness), projected central dispersion, 
and central surface brightness (from eq. [17] of Paper I); ro and fg are the King core radius ro = [9<7 2 /(47rGp c )] 1//2 and the fraction of the total stellar mass 
inside ro; T is the apparent central mass-to-light ratio [see eq. Ell: bidm and r;, DM are mass and half-mass radius of gravitationally bound DM; r p is the radius 
where density of DM and stars becomes equal; r,„ is the radius where enclosed mass of DM becomes equal to that of stars. All parameters are averaged over 
t = 2.67 ... 3 Gyr. There is no data for the model DO;, as in this model the subhalo is completely disrupted by t ~ 0.79 Gyr. 
"Parameters for the moment of time t ~ 0.28 Gyr. 



ISpitzerl l ll987l) . The most sensitive dynamic age para meter is 
the core density p c , which according to the model of Spitzer 
lll987l) evolves with time t as p c oc (1 -f/f C oii)" I//() ' 86 , where Wi 
is the core collapse time. For our model S, p c (3 Gyr)/ j o e (0) = 
1300/360 (see Table 0, so t/t coU ~ 2/3 at t = 3 Gyr. Galac- 
tic GCs are known to span the whole spectrum of dynamic 
ages, ranging from dynamically young with f/f co u <C 1 (such 
as to Centauri and NGC 2419) to post-core-collapse systems 
with t/t coU ^ 1 (~ 20% of all Galactic GCs). As you can 
see, at t $C 3 Gyr our models span a large range of dynamic 
ages, corresponding to many (probably most) Galactic GCs. 
In addition, the ratio of the predicted model core collapse time 
f co n ~ 4.5 Gyr to the model orbital time f or b ~ 0.3 Gyr is ap- 
proximately the same as for dynamically old Galactic GCs, 
which have f co n/f or b ~ 15 (assuming that f co u ~ 12 Gyr and 
forb ~ 0.8 Gyr). As a result, our models can correctly describe 
both secular evolution and tidal stripping of many Galactic 
GCs. 

One could treat the observed long-term dynamic evolution 
of stars in our models as a numeric artifact and a nuisance. 
Instead, by arguing that this effect reflects main features of 
gravothermal instability in real GCs, we can use it to explore 
the impact of long-term dynamic evolution on the properties 
of hybrid (stars + DM) GCs. 

3.2. Isolated Models 

In Paper I we demonstrated that a warm collapse of a ho- 
mogeneous isothermal stellar sphere inside a live DM halo 
(either NFW or Burkert) produces a GC-like cluster with an 
outer density cutoff in the stellar density distribution resem- 
bling a tidal feature in King model. Our results confirmed the 
idea of Peebles ( 1984) that a presence of DM can be an alter- 
native explanation for apparent "tidal" radial density cutoffs 
observed in some GCs. Here we use models S and SD„£ to 
check if this explanation can be extended to stellar clusters 
which have experienced significant secular evolution. 

In Figure |2] we show the stellar radial density profiles for 
models S (short-dashed lines) and SD„ (solid lines) for the 
same four time intervals as in Figure^ As you can see, flat- 
tening of the outer density profile caused by the dynamic evo- 



lution of the cluster does not remove an apparent cutoff fea- 
ture observed in DM-dominated GCs. As time goes on, the 
slope of the outer density profile becomes more shallow for 
both purely stellar and hybrid GCs, but the relative change of 
the slope caused by the presence of DM stays approximately 
the same. Also, the radius where the two profiles start di- 
verging, stays approximately the same (~ 18 pc) for the NFW 
profile, and gradually increases (from ~ 35 to ~ 45 pc) for 
the Burkert profile. For both types of DM halos, this radius is 
very close to the radius r,„ where the inclosed masses of stars 
and DM become equal (see Table|2j. 

The reason for the persistence of the density cutoff in the 
course of the secular evolution of our cluster is the same as 
for the appearance of the cutoff at the end of the initial violent 
relaxation phase (see Paper I). The original cutoff is caused 
by the fact that at large radii (r > r,„) the potential of the hy- 
brid GC is dominated by DM. As a result, a smaller fraction 
of stars, ejected from the cluster during the violent relaxation, 
can populate outer halo, creating a density cutoff at r ~ r m . 
Similarly, during long-term dynamic evolution of the cluster, 
caused by encounters between individual stars, stars are be- 
ing ejected from the core (core evaporation), with very few of 
them populating the halo beyond r m because of the potential 
being dominated by DM at such large radii. 

We use isolated DM-only models D„ to see how numeri- 
cal artifacts affect the DM density distribution in our simula- 
tions. lHavashi et alJ (120031) discussed the impact of discrete- 
ness effects of A^-body simulations on radial density profiles 
of isotropic NFW halos. According to these authors, for the 
first ^100 crossing times at the virial radius, the density at the 
center of the halo is decreasing due to heating by faster mov- 
ing particles, causing the core to expand. At the end of this 
phase, the central velocity dispersion becomes comparable to 
the velocity dispersion at the scale radius r s . (In NFW mod- 
els the velocity dispersion is highest around r s , and becomes 
smaller for both smaller and larger radii.) 

We see a similar trend in our NFW model, D„. In our case, 
the total evolution time is ~ 24 crossing times at the virial ra- 
dius T v j r =124 Myr, so at the end of the simulations the model 
is still in the core-heating regime. As you can see in Figure^, 
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FIG. 2.— Averaged stellar radial density profiles for four time intervals: (a) t = 0.12... 0.37 Gyr; (b) t = 0.37... 0.87 Gyr; (c) t = 0.87... 1.61 Gyr; (d) 
t = 1.61 ... 3 Gyr (the same time intervals as in Figure^. In the bottom-right parts of the panels, the sequence of models is as follows (from left to right): SD„ 
(solid line; colored green in electronic edition), SDO„ (dotted line; colored green in electronic edition), SD/, (solid line; colored blue in electronic edition), SDOj 
(dotted line; colored blue in electronic edition), and S (short-dashed line). 



at t ~ 3 Gyr the inner DM density profile is reasonably close 
to the initial one down to a radius of ~ 2.5 pc ~ 1 .7£dm- 

In the case of the Burkert profile, the initial central velocity 
dispersion is also smaller than the dispersion at r s , though the 
difference is much smaller than for the case of NFW halo. In 
our model D/,, we see a slight decrease in the DM density in 
the inner part of the halo at the end of the simulations (Fig- 
ure|2j>), with a reasonably accurate profile down to a radius of 
-3.5 pc~2.3e D M- 

NFW and Burkert models would be in equilibrium only if 
they had infinitive size and mass. As our models are truncated 
at at a finite radius r v ; r , the DM density in the outer parts of the 
halos becomes lower with time. As you can see in Figure|3j 
even at the end of the simulations the deviation of the outer 
DM density profile from the theoretical one is not significant 
within r v ; r , and is negligible for r < 3r s . Most of the changes 
in the outer density profiles happen within first one or two T YiI , 



which is comparable to the orbital period in models DO„,6 and 
SDO„> All the above let us conclude that the truncation of 
our models at the radius of r wlI should not have a noticeable 
impact on the results of our simulations. 

3.3. Evolution in Tidal Field: Preliminary Analysis 

lHavashi etafl J2001 showed that an isotropic NFW sub- 
halo, following a circular orbit inside the static potential of a 
host galaxy with NFW profile, can be completely disrupted 
by tidal forces after a few orbits if its tidal radius r^ is 
smaller than — 2r s . These authors argued that the relative 
easiness to disrupt an NFW satellite is linked to the fact that 
when an NFW halo is instantaneously truncated at a radius 
f < r bind — 0.11 r s , it becomes unbound (with the total energy 
E tot of the remnant becoming positive). The situation with 
singular isothermal spheres is completely different, as such 
halos stay gravitationally bound for any r, and potentially can 
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FIG. 3. — Averaged DM radial density profiles for the time interval t = 1.61 ... 3 Gyr (the same time interval as in Fig. 01). (a) NFW models, (b) Burkert 
models. Long-dashed lines show analytical profiles for undisturbed DM halos. Thick solid lines correspond to models D„ j, (no stars, no static potential), thin 
solid lines show the profiles for models SD„ 4 (DM + stars, no static potential), and short-dashed lines correspond to models SDO„ j, (DM + stars on orbit in static 
potential). Vertical long-dashed, dotted, and solid lines mark the values of eoM. '".1, and r v ; r , respectively. 




FIG. 4. — Total energy E tot of halos instantaneously truncated at a ra- 
dius r. Here r is in r s units, and £ t ot is in GmJ ir /r, units. Thick and thin 
lines (colored green and blue in electronic edition, respectively) correspond 
to NFW and Burkert profiles, respectively. Long-dashed lines correspond to 
the initial configuration of DM-only D„j, models, and solid lines correspond 
to sufficiently relaxed (f = 0.25 Gyr) models SD„ j, which have a stellar core. 



FIG. 5. — Radial profile of tidal acceleration dax/dR for NFW (thick solid 
line; colored green in electronic edition) and Burkert (thin solid line; colored 
blue in electronic edition) halos with a concentration of C = 6.75. Here R is 
in R s units, and da^/dR is in GM yiI /R^ units. Vertical dotted lines show the 
pericentric and apocentric distances for the subhalo orbit in models DO„ /, 
and SDO,, ;,. Horizontal short-dashed line corresponds to zero radial tidal 
acceleration. 



survive indefinitely in an external tidal field. 

In Figure|4]we compare the binding properties of our NFW 
and Burkert halos (thick and thin long-dashed lines, respec- 
tively). As you can see, Burkert halos are similar to NFW 
ones in becoming unbound if truncated below a certain ra- 
dius. Burkert halos appear to be much easier to disrupt tidally 
than NFW halos: in the case of Burkert profile, rbind — 1 -66r s , 
which is ~ 2.1 times larger than for NFW profile. Also, the 
positive total energy of a Burkert halo truncated to r < rt,md 
is ^3.3 times larger than the corresponding quantity for an 



NFW halo with the same concentration and mass (see Fig- 
ure|4j. 

In Figure 13 we compare the strength of the tidal force as 
a function of radius for two types of host galaxies: with 
NFW profile (thick line) and with Burkert profile (thin line). 
More specifically, in this figure we plot the radial dependence 
of the radial gradient of gravitational acceleration da^/dR = 
d/dR[GM(R)/R 2 ]. A product of this quantity on the linear 
size of a subhalo gives an estimate of the differential (tidal) 
acceleration between two opposite parts of the subhalo. Two 
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FIG. 6. — Time evolution of the gravitationally bound DM mass of subhalos in models SDO„ j, (thick lines) and DO„j, (thin lines), (a) NFW halos (models 
SDO„ and DO,,), (b) Burkert halos (models SDOj and DO/,). Short-dashed lines show the evolution of the bound stellar mass in the corresponding SDO models. 
Dotted line in the panel (b) shows the evolution of DM bound mass for subhalo from the model SDO4 which at t = 2.13 Gyr was removed from the static potential 
of the host galaxy and was let to evolve in isolation. 



vertical dotted lines in Figure |5] show the range of radial dis- 
tances covered by our models. 

As you can see in Figure [5] for NFW profile the quan- 
tity da/i/dR is always positive (meaning that the tidal force 
is always stretching a subhalo in radial direction). In the 
small radii limit, da^jdR asymptotically approaches a con- 
stant, which is equal to 2/3 [ln(l +C)-C/(1 +C)]" 1 GM^jR\. 
For Burkert halos, the radial tidal force reaches a maximum at 
R = 1 .76328/?,, becomes equal to zero at R = 0.96340fl s , and 
asymptotically approaches a negative constant -2/3[ln(l + 
Q+l/2 ln( 1 + C 2 ) - arctanC]" 1 GM^/R\ at small radii. The 
negative sign for da^/dR means that at R < R s the radial tidal 
force becomes compressing instead of being stretching. In the 
interval of radii covered by the orbital motion of our subha- 
los, the NFW host halo has significantly stronger radial tidal 
acceleration than the Burkert halo (see Figure[f>J. 

To summarize the above analysis, the Burkert halos are eas- 
ier to disrupt tidally than NFW halos. On the other hand, 
a tidal field of a Burkert host galaxy appears to be less dis- 
rupting than that of an NFW host galaxy with the same mass 
and concentration. In addition, Burkert halos have an unusual 
property of having a compressing radial tidal force within the 
scale radius R s . One has thus to resort to numerical simula- 
tions to understand differences in substructure evolution for 
NFW and Burkert cases. 

3.4. Tidal Stripping of DM-ordy subhalos 

In our models DO„i, a DM-only (NFW or Burkert) subhalo 
is orbiting on eccentric (R a /R p = 5) orbit inside a static poten- 
tial of the host (N FW or Burkert) gala xy. We use the same 
definition of r U d as Havas hTet al.N2003T) : 



m(hid) 



' tid 



R dM] M(R) 



M(R) dR J R 



(1) 



Here M(R) and m(r) is enclosed mass as a function of radius 
for the host and satellite halos, respectively. At R = R p , the 
subhalos have the following values of tidal radii: r U d = 314 pc 



= 1.73r. s for NFW profile, and r tid = 498 pc = 2.15r s for Burk- 
ert profile. 

In Figure [6^, we show the evolution of the gravitationally 
bound DM mass for model DO,, (solid thin line). As you 
can see, tidal stripping is severe in this model, with only 
~ 3.8% of the total mass surviving as a bound structure after 
11 orbits (see Table 0. This behavior is similar to the criti- 
cal case of Havashi et al. (2003) with r t a/r s = 2.1 (see their 
Fig. 7), which separates their models staying relatively intact 
(with r t id//".s > 2.1) and completely disrupted models (with 
ftid/r s < 2.1). In our case, this ratio is slightly smaller than 
2. 1 : r t id/r s = 1 .73. Slightly larger resilience to tidal disruption 
exhibited by o ur mod el DO„ can be explained by the fact that 
lHavashi et alJ (12003) used "local Maxwellian approximation" 
to set up the initial particle distribution in their models, which 
according to Kazantzidis et al. (2004) leads to artificially high 
disruption rate for NFW subhalos. 

Even for the very last snapshot of model DO,, (at t = 3 Gyr), 
the stripped-down subhalo appears to be reasonably stable 
against total disruption by the tidal forces. Indeed, the bind- 
ing radius of the remnant rb; n d = 15 pc is significantly smaller 
than the current value of the tidal radius r t ;d = 56 pc even at 
the very end of the simulations. (To estimate r^, we use the 
actual value of the pericentric distance of R p = 413 pc which 
is smaller than the original value of R p = 600 pc because of 
the dynamic friction experienced by the subhalo moving in 
the halo of unbound DM particles.) We expect thus that our 
NFW subhalo should survive as a gravitationally bound rem- 
nant for a few more orbits (and perhaps indefinitely). 

The fate of a Burkert DM subhalo orbiting inside a Burkert 
host galaxy (model DOb) is completely different, as can be 
seen in Figure [6j> (solid thin line). After only 3 pericentric 
passages, the subhalo ceases to exist as a coherent structure. 
To make sure that the failure to find a bound structure after 
t = 0.79 Gyr is not an artifact of our halo finding algorithm, 
we measure both binding radius rb; n d = 625 pc and tidal radius 
r t id = 175 pc for the last snapshot containing a bound structure. 
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FIG. 7. — Grey-scale maps of DM column density for models SDO„ (two left panels) and SDO4 (two right panels) for two moments of time: t ~ 0.43 Gyr 
(when the subhalos are approaching the pericenter for the second time; two top panels) and t = 3 Gyr (at the end of simulations; two bottom panels). The orbit of 
the subhalo is seen face-on. In all panels we use the same logarithmic scale for the values of DM column density. Cross marks the location of the center of the 
host galaxy, and two concentric circles correspond to apocentric and pericentric distances of the subhalo's orbit. 



As you can see, at this point the subhalo is bound to be dis- 
persed by tidal forces as its current tidal radius is significantly 
(almost 4 times) smaller than the binding radius. 

Two main conclusions we arrive at in this section are (1) 
our NFW DM-only subhalo is relatively resilient to total tidal 
disruption, and (2) a Burkert DM subhalo orbiting inside a 
Burkert host galaxy is much easier to disrupt that for the case 
of an NFW DM subhalo orbiting inside an NFW host. 

3.5. Tidal Stripping of Hybrid GCs 

We are now turning to the results of our two principle mod- 
els: SDO„ and SDO/,. In these models, a hybrid GC (stars 
+ DM) is orbiting inside the same host galaxies and along 
the same eccentric orbit as in models DO,,./,. Before being 
placed in the static potential of the host galaxy, we allow stars 
and DM to relax in isolation for 120 Myr (model SDO„) and 
170 Myr (model SDO,,). 



As you can see in Figure[6] presence of a stellar core, with a 
mass of mere ~ 0.9% of the total mass, inside a DM subhalo 
makes the subhalo much more resilient to tidal forces disrup- 
tion. In the case of NFW profile (Figure^; see also Tabled, 
the DM mass of the tidally stripped remnant at t ~ 3 Gyr is 
7 times larger in the presence of the stellar core than in the 
DM-only case. It is also 3 times larger than the mass of the 
stellar core. For Burkert profile, the difference is even more 
dramatic: whereas a starless Burkert subhalo is completely 
disrupted after 3 orbits, the subhalo with a stellar core sur- 
vives till the end of the simulations, with ~ 32% of the total 
mass of the remnant being in DM form. 

Figure ^provides an exp lana tion for such a marked differ- 
ence. As we discussed in § 13.31 both NFW and Burkert halos 
become unbound if instantaneously truncated below a certain 
radius (0.77r s for NFW profile, and \.bbr s for Burkert pro- 
file). As you can see in Figure|4] the presence of a stellar GC- 
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FIG. 8. — Final radial density profiles for (a) NFW models and (b) Burkert models (for the time interval ; = 2.67 ... 3 Gyr). Thick and thin (colored red in 
electronic edition) lines correspond to DM and stars, respectively. Long-dashed lines show the profiles of DM in the absence of stars and tidal field (analytical 
models) and stars in the absence of DM (model S). Solid lines correspond to SDO„ /, models. Vertical long-dashed, short-dashed, dash-dotted, and dotted lines 
mark the values of edm, T m , current tidal radius for SDO models rga, and r s . 




FIG. 9. — Final radial profiles of observable quantities (for the time interval ( = 2.67 ... 3 Gyr - the same as in Fig.|S|. (a) V-band surface brightness Ey. (b) 
Line-of-sight stellar velocity dispersion o>. Thick and thin solid lines (colored green and blue in electronic edition, respectively) correspond to models SDO„ 
and SDO/,, respectively. (There are two DSO/, profiles: the upper one corresponds to stars along the major axis of the cluster, and the lower one is for stars in 
the plane perpendicular to the major axis.) Long-dashed lines show the profiles for the model S. Thick vertical short-dashed line marks the value of r m for model 
SDO,, (see Tablel2l. Vertical dash-dotted lines correspond to the current tidal radius for SDO models r t a. To make these plots we use the projection technique 
described in Appendix B of Paper I. 



like cluster at the center of either halo makes it significantly 
more bound. In the case of NFW profile (thick solid line), 
the halo becomes bound for virtually any truncation radius r 
(being similar in this regard to a singular isothermal sphere). 
For Burkert halo (thin solid line), there is still a range in trun- 
cation radius r where the halo is formally unbound, but the 
maximum positive value of E tot is significantly lower in the 
presence of a stellar core. Moreover, for r < 0.32r s the halo 
again becomes bound (see Figure^). 

After 7 orbits or ~ 2 Gyr, the mass of DM gravitation- 
ally bound to the remnant in model SDO/, becomes smaller 



than the mass of the stellar core. At this point, the rate of 
the mass loss due to tidal stripping becomes smaller, but still 
non-negligible (see Figure|U>). We tested the possibility that 
the observed decrease in bound DM mass after t ~ 2 Gyr is a 
numerical artifact, caused by a small number (< 10 4 ) of DM 
particles attached to the remnant, by resimulating the late evo- 
lution of the subhalo in the absence of the static gravitational 
potential of the host galaxy. For this, we used the gravitation- 
ally bound remnant from the t = 2.13 Gyr snapshot of model 
SDO/,. The bound DM mass evolution for this additional run 
is shown with dotted line in Figure |6ji. As you can see, in the 
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absence of external tidal field the bound DM mass of the rem- 
nant stays virtually constant. It appears thus very unlikely that 
the late time DM stripping observed in model SDO/, is caused 
by numerical artifacts (such as evaporation of DM particles 
due to two-body interactions). 

In Figure we show DM column density maps for models 
SDO,, (left panels) and SDO/, (right panels), for two moments 
of time: after ~ 1 .5 orbits (top panels) and after 10-11 orbits 
(bottom panels). Two concentric circles show the range of ra- 
dial distances covered by the models. At the earlier moment 
of time, you can see a long stream of DM particles, stripped 
during the first pericentric passage, with a relatively massive 
bound structure in the center of the stream. At the end of 
the simulations, multiple streams of stripped particles mix to- 
gether to create a fuzzy thick disk of DM. The gravitationally 
bound remnant is barely visible at this point (especially for 
Burkert model). 

In Figure [3] we show averaged DM radial density profiles 
for t = 1 .61 ... 3 Gyr for models D,,./, (thick solid lines), SD„ /, 
(thin solid lines), and SDO,,,/, (short-dashed lines). As you 
can see, the DM density in the innermost part of the NFW 
halo is not modified by external tidal field. Outside of the 
stars-dominated central region, the DM density profile be- 
comes significantly steeper in the presence of tidal field. At 
this point, in models SDO,,./, DM is stripped down to the orig- 
inal scale radius r s . 

Analysis of Table|2]shows that the parameters of our stellar 
clusters at the end of the simulations are very similar for all 
our models containing stars (S, SD„,/,, and SDO,,./,). For SD„ 
and SDO„ models, the presence of relatively large amounts 
of DM at the center of the cluster leads to somewhat larger 
values of the central dispersion a c and smaller values of the 
King core radius ro. The value of the apparent mass-to-light 
ratio T, defined as (see Paper I) 



T = 



9oqTgc 
2irG( R 



(2) 



hb 



where Tec — 1 -45 is the assumed baryonic mass-to-light ratio 
in GCs and £o is the projected stellar surface mass density at 
the center of the cluster, is larger by ~ 18% for a stellar cluster 
in NFW halo. Overall, the presence of an external tidal field 
does not seem to change the global structural and dynamic 
parameters of GCs with DM in a noticeable way. 

As can be seen in Figure|2] tidal stripping does not modify 
significantly the stellar radial density profiles of hybrid GCs: 
at any moment of time, density profiles for SDO,,./, models 
(dotted lines) are very close to the density profiles of SD,,,/, 
models (solid lines). 

In Figure[8]we show the final radial density profiles for both 
stars and DM for models SDO„ and SDO/,. As you can see, 
DM still dominates stars in the outskirts of the stellar cluster. 
In the stars-dominated area, the DM density profile is steeper 
because of the adiabatic contraction of DM in the presence of 
stars. DM density profiles are dramatically modified both by 
the presence of a stellar core and by tidal stripping. 

In model SDO/,, stellar cluster becomes tidally limited by 
the end of the simulations. By t = 3 Gyr, around 0.6% of 
stars (or ~ 60 stellar particles) are not gravitationally bound 
to the cluster, forming distinctive trailing and leading stellar 
tidal tails. (Conversely, not a single stellar particle has been 
tidally stripped by t = 3 Gyr in model SDO„.) In model SDO/,, 
the shape of the cluster becomes increasingly non-spherical at 
large radii at t ~ 3 Gyr. As can be seen in Figure |9^, there 



is an excess of stars in the outskirts of the cluster along its 
major axis, and a sharp tidal cutoff around the analytical tidal 
radius near the plane perpendicular to the major axis. (To 
calculate surface brightness and velocity dispersion profiles 
for SDO/, model, we used all stellar particles - both bound 
and unbound.) 

As Figure^ shows, the final surface brightness profiles of 
the stellar clusters in SDO,, /, models look remarkably simi- 
lar to the corresponding profiles of Galactic GCs. (We as- 
sumed that the baryonic V-band mass-to-light ratio in GCs is 
T"gc = 1.45 - see discussion in Paper I.) The "dent" feature 
observed in the surface brightness profile of a freshly relaxed 
hybrid (DM + stars) GC (see Paper I) has been removed in 
mode ls SDO,,./, by secular evolution of the stellar cluster (see 
§ 13. 1> . The apparent "tidal" cutoff in the outer surface bright- 
ness profile of model SDO„ is caused by the presence of sig- 
nificant amounts of DM in the outskirts of the cluster (simi- 
larly to Paper I). In the case of Burkert halo (model SDO/,), 
the cutoff is of truly tidal nature. 

In Figure^ we show the final line-of-sight velocity disper- 
sion profiles for the stellar clusters in SDO,,,/, models. In the 
case of NFW halo (thick solid line), the line-of-sight velocity 
dispersion appears to be uniformly inflated by a small fac- 
tor (~ 10%) across all radii, which can be misinterpreted as a 
purely stellar cluster with a somewhat larger value of baryonic 
mass-to-light ratio. For the Burkert halo (thin solid lines), the 
radial velocity dispersion profile for the stars near the plane 
perpendicular to the major axis is very close to the profile of 
a purely stellar cluster, whereas the stars along the major axis 
show signs of being tidally heated. 

4. CONCLUSIONS 

DM subhalos with Burkert density profile are much easier 
to disrupt tidally in the potential of the host galaxy than sub- 
halos with NFW profile. We link this effect to the difference 
in the binding properties of both types of halos, with the to- 
tal energy of the Burkert halos becoming positive if truncated 
at a significantly larger radius than for NFW halos. Setting 
a low-mass (~ 1% of the total mass) dense stellar core at the 
center of either NFW or Burkert halo makes them much more 
resilient to tidal disruption. 

Primordial GCs with NFW DM halo can survive severe 
tidal stripping in the host galaxy, with DM still being the dom- 
inant mass component (though not by a large margin) in the 
tidally stripped-down remnant. DM is concentrated in the out- 
skirts of the remnant. As a result, an apparent "tidal" cutoff 
in the radial surface brightness profile in isolated warm col- 
lapse models (Paper I), caused by the presence of DM, is also 
present in our tidally stripped NFW model. 

We used warm collapse hybrid models to show that neither 
secular evolution of the stellar cluster nor severe tidal strip- 
ping change noticeably the inferred core mass-to-light ratio 
T. For both flat-core and cuspy DM halo profiles, T stays 
close to the purely baryonic value. 

Tidal stripping can remove almost all DM from primordial 
GCs with Burkert DM halo. The remaining DM is dynam- 
ically unimportant, and cannot prevent stars being stripped 
off by tidal forces of the host galaxy. This result makes GCs 
possessing obvious tidal tails (the only known example being 
Palomar 5) be fully consistent with primordial scenarios of 
GC formation. 

Secular evolution of a DM-dominated GC does not change 
the main results of Paper I derived for freshly relaxed warm- 
collapse systems. In particular, in both unevolved and evolved 
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clusters, presence of significant amounts of DM in the out- 
skirts of the cluster manifest itself as a "tidal" cutoff in the 
outer part of the radial surface brightness profile. It is not 
clear if the "extratidal" features of cold collapse models from 
Paper I will be pr eserv ed in significantly evolved clusters. As 
we discussed in § 12.11 simulating a cold-collapse hybrid GC 
to address this issue is not feasible with the present day tech- 
nology. 

The above results reinforce our conclusion from Paper I that 
a presence of obvious tidal tails is probably the only obser- 
vational evidence which can reliably rule out a presence of 
significant amounts of DM in GCs. (But it cannot rule out 
primordial scenarios of GC formation.) 

To summarize, the results presented in both Paper I and this 



paper suggest that the whole range of features seen in Milky 
Way GCs (from apparently truncated profiles of some clusters 
to the extended tidal tails of Palomar 5) can be consistent with 
the primordial picture of GC formation, given that there was a 
range of the initial virial ratios for stellar clusters (from "cold" 
to "hot" collapses), tidal stripping histories, and/or inner DM 
density profiles. 
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